module solvers
   !=============================================================================
   ! Aims:
   !   1) Given the coefficients, sources and a initial guess, solve
   !      with the msi method the system of algebraic equations like AX = B,
   !      where A has 5 non-null diagonals for pl and
   !                  9 non-null diagonals for u, v and T
   !   2) Calculate the L1 norm of the residual of the algebraic equations,
   !      for 5 and 9 non-null diagonals, considering only real volumes.
   !
   ! Latest version: 14 Oct 2012
   !
   ! Subroutines given by this module:
   !   1) norma_L1_5d
   !   2) norma_L1_9d
   !
   ! Subroutines given by module msi2d5
   !   3) lu2d5
   !   4) fb2d5
   !
   ! Subroutines given by module msi2d9
   !   5) lu2d9
   !   6) fb2d9
   !=============================================================================

   use msi2d5
   use msi2d9

   implicit none

contains

   !=============================================================================

   subroutine norm_l1_5d(nx, ny, var, b, a, norm)
      implicit none
      integer, intent(in) :: nx   ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny   ! Number of volumes in eta direction (real + fictitious)
      real(8), dimension(nx*ny), intent(in) :: var ! Unknown variable
      real(8), dimension(nx*ny), intent(in) :: b   ! Source vector of the linear system
      real(8), dimension(nx*ny,5), intent(in) :: a   ! Coefficients of the linear system
      real(8), intent(inout) :: norm

      ! Auxiliary variables
      integer :: i, j, np, nps, npn, npw, npe, aux_i

      ! Norm is calculated taking into account only real volumes
      do j = 2, ny-1
         aux_i = nx*(j - 1)
         do i = 2, nx-1
            np = aux_i + i
            nps = np - nx
            npn = np + nx
            npw = np - 1
            npe = np + 1
            norm = norm + dabs(a(np,1)*var(nps) + a(np,2)*var(npw) &
               + a(np,3)*var(np) + a(np,4)*var(npe) + a(np,5)*var(npn) - b(np))
         end do
      end do

   end subroutine norm_l1_5d

   !=============================================================================

   subroutine norm_l1_9d(nx, ny, var, b, a, norm)
      implicit none
      integer, intent(in) :: nx   ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny   ! Number of volumes in eta direction (real + fictitious)
      real(8), dimension(nx*ny), intent(in) :: var ! Unknown variable
      real(8), dimension(nx*ny), intent(in) :: b   ! Source vector of the linear system
      real(8), dimension(nx*ny,9), intent(in) :: a   ! Coefficients of the linear system
      real(8), intent(inout) :: norm

      ! Auxiliary variables
      integer :: i, j, np, nps, npn, npw, npe, npsw, npse, npnw, npne, aux_i

      ! Norm is calculated taking into account only real volumes

      do j = 2, ny-1
         aux_i = nx*(j - 1)
         do i = 2, nx-1

            np = aux_i + i
            nps = np - nx
            npn = np + nx
            npw = np - 1
            npe = np + 1
            npsw = nps - 1
            npse = nps + 1
            npnw = npn - 1
            npne = npn + 1

            norm = norm + dabs(a(np,1)*var(npsw) + a(np,2)*var(nps) &
               + a(np,3)*var(npse) + a(np,4)*var(npw) + a(np,5)*var(np) &
               + a(np,6)*var(npe) + a(np,7)*var(npnw) + a(np,8)*var(npn) &
               + a(np,9)*var(npne) - b(np))

         end do
      end do

   end subroutine norm_l1_9d

   !============================================================================

   subroutine tdma(a, b, x, n)
      ! O TDMA resolve sistemas tridiagonais do tipo Ax = b de ordem n
      implicit none
      integer n, i
      real*8, dimension(:) :: b(n), x(n), p(n), q(n)
      real*8, dimension(:,:) :: a(n, n)

      p(1) = -a(1,2)/a(1,1)
      q(1) = b(1)/a(1,1)

      do i = 2, n
         if (i < n) p(i) = -a(i, i + 1)/(a(i, i) + a(i, i-1)*p(i-1))
         q(i) = (b(i) - a(i, i-1)*q(i-1))/(a(i, i) + a(i, i-1)*p(i-1))
      end do

      x(n) = q(n)

      do i = n-1, 1, -1
         x(i) = p(i)*x(i + 1) + q(i)
      end do

   end subroutine tdma

   !============================================================================

   !> Solves a tri-diagonal linear system
   subroutine tdma3d(n, a, b, x)
      implicit none
      integer, intent(in) :: n      !< Number unknowns
      real(8), intent(in) :: a(n,3) !< Tri-diagonal matrix
      real(8), intent(in) :: b(n)   !< Source
      real(8), intent(out) :: x(n)  !< Solution

      ! Auxiliary variables
      integer :: i
      real(8), dimension(n) :: P
      real(8), dimension(n) :: Q

      i = 1

      P(i) = -a(i,3)/a(i,2)

      Q(i) = b(i)/a(i,2)

      do i = 2, n

         P(i) = -a(i,3)/(a(i,2) + a(i,1)*P(i-1))

         Q(i) = (b(i) - a(i,1)*Q(i-1))/(a(i,2) + a(i,1) * P(i-1))

      end do

      i = n

      x(i) = Q(i)

      do i = n-1, 1, -1

         x(i) = x(i+1)*P(i) + Q(i)

      end do

   end subroutine tdma3d

   !============================================================================

   !> \brief Calculates the coefficient for convergence criteria of a 5-diagonal matrix.
   !! This coefficient is calculated only for real volumes.
   real(8) function get_maxc_5d(nx, ny, a)
      implicit none
      integer, intent(in) :: nx   !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny   !< Number of volumes in the eta direction (real+fictitious)
      real(8), dimension (nx*ny,5), intent(in) :: a   !< Coefficients of the linear system

      ! Inner variables
      integer :: i, j, np
      real(8) :: c, maxc

      ! Looking for the maxval of c = ( sum_nb | A_nb | ) / | A_P |

      maxc = 0.d0
      do i = 2, nx-1
         do j = 2, ny-1
            np   = nx * (j-1) + i
            c = (sum(abs(a(np,1:2))) + sum(abs(a(np,4:5))))/abs(a(np,3))
            if (maxc < c) maxc = c
         end do
      end do

      get_maxc_5d = maxc

   end function get_maxc_5d

   !============================================================================

   !> \brief Calculates the coefficient for convergence criteria of a 9-diagonal matrix.
   !! This coefficient is calculated only for real volumes.
   real(8) function get_maxc_9d(nx, ny, a)
      implicit none
      integer, intent(in) :: nx   !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny   !< Number of volumes in the eta direction (real+fictitious)
      real(8), dimension(nx*ny,9), intent(in) :: a   !< Coefficients of the linear system

      ! Inner variables
      integer :: i, j, np
      real(8) :: c, maxc

      ! Looking for the maxval of c = (sum_nb|A_nb|)/|A_P|

      maxc = 0.d0
      do i = 2, nx-1
         do j = 2, ny-1
            np   = nx*(j - 1) + i
            c = (sum(abs(a(np,1:4))) + sum(abs(a(np,6:9))))/abs(a(np,5))
            if(maxc < c) maxc = c
         end do
      end do

      get_maxc_9d = maxc

   end function get_maxc_9d

   !============================================================================

end module solvers
